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Higher-order perturbative calculations in Quantum (Field) Theory suffer from 
the factorial increase of the number of individual diagrams. Here I describe an 
approach which evaluates the total contribution numerically for finite temper- 
ature from the cumulant expansion of the corresponding observable followed 
by an extrapolation to zero temperature. This method (originally proposed by 
Bogolyubov and Plechko) is applied to the calculation of higher-order terms 
for the ground-state energy of the polaron. Using state-of-the-art multidimen- 
sional integration routines 2 new coefficients are obtained corresponding to a 
4- and 5-loop calculation. 
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1. Introduction 

Highly accurate measurements require precise theoretical calculations which 
perturbation theory can yield if the coupling constant is small. However, 
in Quantum Field Theory (QFT) the number of diagrams grows factorially 
with the order of perturbation theory and they become more and more com- 
plicated. The prime example is the anomalous magnetic moment of the elec- 
tron where new experiments^ need high-order quantum-electrodynamical 
calculations but the number of diagrams for them "explodes" as shown by 
the generating function 

r(a) = l + a + 7a 2 + 72a 3 + 891a 4 + 12672a 5 + 202770a 6 + ... (1) 

There are ongoing efforts^ to calculate all 12672 diagrams in 0(a 5 ) - a 
huge, heroic effort considering the complexity of individual diagrams and 
the large cancellations among them. 

Obviously new and more efficient methods would be most welcome for 
a cross-check or further progress. 
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2. A new method (applied to the polaron g.s. energy) 

Here I present a "new" method which - as I learned during the conference 
- was already proposed 20 years by Bogolyubov (Jr.) and Plechko (BP) ^. 
However, to my knowledge it has been never applied numerically which 
turned out to be quite a challenging task. 

The BP method is formulated for the polaron problem, a non-relativistic 
(but non-trivial) field theory describing an electron slowly moving through 
a polarizablc crystal. Due to medium effects its energy is changed and it 
acquires an effective mass : E p = Eq + p 2 /(2m*) + . . .. The aim is to 
calculate the power series expansion for the g.s. energy 

E (a) =: ^e„ a " (2) 

71=1 

as function of the dimensionless electron-phonon coupling constant a. The 
lowest-order coefficients are well-known^ ( e\ — — 1 , = —0.01591962 ) 
but since Smondyrev's calculation" in 1986 

e 3 = -0.00080607 (3) 

there has been no progress towards higher-order terms. 

This will be remedied by the first numerical application of the BP 

method. For this purpose the path integral formulation of the polaron prob- 

7 

lcm will be used where the phonons have been integrated out exactly . For 
large Euclidean times (3 this gives the following effective action 

rl rl f ,l. 7 a f 13 , /■* tM fd 3 k explik- <x(t) - x(i'))] 

(4) 

which will be split into a free part So and an interaction term Si. The g.s. 
energy may be obtained from the partition function 

Z(f3) = jv^xe- 8 ^ e-P E ° (5) 

at asymptotic values of (3, i.e. zero temperature. The central idea is to use 
the cumulant expansion of the partition function 



Z{f3) = Zo exp 



n—1 



(6) 



where the A„(/3)'s are the cumulants w.r.t. 5i. These are obtained from the 
moments 

m n := Af £v 3 x (S 1 [x] ) n e" 5 "^ 1 , m = 1 (7) 
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by the recursion relation (see, e.g. Eq. (51) in Ref. 8) 



A n +i — m n+ i 



n-l 



k=o ^ k 



(8) 



Explicitly the first cumulants read 

Ai = mi , A 2 = m 2 — m\ , A3 = m 3 — 3 m 2 mi + 2 m\ 

A4 = m4 — 4 m3 mi — 3 m 2 + 12 to 2 m 2 — 6 m^ (9) 

A5 = m 5 — 5 m4 mi — 10 m 3 m 2 + 20 m 3 m 2 + 30 m 2 mi — 60 m 2 m^ + 24 mf 

By construction m„ cx a™ and Eq. (8) shows that the cumulants share this 
property. Thus we immediately obtain 



1 f-)™ +1 



(10) 



The functional integral for the moments can be done since it is Gaussian 
The integrals over the phonon momenta k m , m = 1 . . . n can also be per- 
formed if the m th propagator is written as 



1 



2 Jo 



du m exp 



--\c 2 v 

~ LL r, 



(11) 



Then one obtains 

_ (-)"a 



(47t)«/ 2 n 

v ' m=l 



dt r 



tm 



dt^ 



du r 



exp 



(im i m ) 



m— 1 



[det^(ti...t n ,t , 1 ...t' n ; «!■■■««)] 3/2 • (12) 



Here the (n x ri)- matrix A 
1 



\ti tj 



\t'-t'A 



m Sij . (13) 



is non-analytic in the times U,^, but analytic in the auxiliary variables Mj. 



3. Numerical procedures and results 

The task is now to perform the (3n)-dimensional integral over ti,t\,Ui for 
large enough [3 in the expression for the cumulants/moments. It is clear 
that any reduction in the dimensionality of the integral will greatly help 
in obtaining reliable numerical results in affordable CPU-time. A closer in- 
spection of the structure of the integrand reveals that 2 integrations over 
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the auxiliary variables (say u n , can always be done analytically. Fur- 

thermore, we do not use Eq. (10) to extract the energy coefficient e n but 



i-r 



a n nl P 



lim 



df3 



= : lim e n (/3) 

p— ►oc 



(14) 



This "kills two birds with one stone" : first the derivative w.r.t. [3 takes away 
one further integration over a time (see Eq. (12) where /3 appears as upper 
limit) requiring that only a (3n — 3)-dimensional integral has to be done 
numerically. Second, it vastly improves the convergence to e„ = e n (fl = oo) 
because now 



-const- 



VP 



-0. 



(15) 



In other words : we obtain an exponential convergence to the value e n 
whereas previously the approach would be very slow, like const//?. This 
exponential convergence of the derivative version has been demonstrated 
analytically for n — 1,2 and numerically for n = 3 (see below). In the 
following we will assume that it holds for all n. After mapping to the hy- 
percube [0, 1] the remaining (3n — 3)-dimensional integral can be evaluated 
by Monte-Carlo techniques utilizing the classic VEGAS program^ or the 
more modern programs from the CUBA library-^. 

We first have tested this approach by comparing with the analytical 
result given in Eq. (3). Fig. 1 shows es(P) and the best fit to the data 
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e 3 =-0.8060(2)10- 3 



Fig. 1. (color online). Monte-Carlo results for the derivative of the 3 rd cumulant as 
function of the Euclidean time j3. The total number of function calls is denoted by ntot 
and the full (open) circles are the points used (not used) in the fit. 
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assuming the /?-dependence (15). Since the asymptotic behaviour is not 
valid for low values of (3 we have eliminated small-/? points successively until 
the resulting x 2 /dof of the fit reaches a minimum. Excellent agreement with 
Smondyrev's result (3) is found. If one allows for a different power of (3 in 
the prefactor of Eq. (15) then the fit gives an exponent —0.55(3) instead 
of —0.5 assumed before. 

However, when extending these calculations to the case n = 4 a very 
slow convergence of the numerical result with the number of function calls 
Utot is observed at fixed f3. Fortunately, a solution was found by perform- 
ing the remaining (n — 2) Mi-integrations not by stochastic (Monte-Carlo) 
methods but by deterministic quadrature rules. This is possible since the 
t^-dependence of the integrand is analytic (see Eq. (13)). We have used 
the very efficient "tanhsinh-integration" method but Gaussian quadra- 
ture is nearly as good. A dramatic improvement in stability results together 
with a reduction of n to t needed for the much smaller values of \e n \ ,n > 3. 
This allows a reliable evaluation of 64 (see Fig. 2 a) and also makes the 
determination of e§ feasible as shown in Fig. 2 b. 



0.54 



0.52 



0.50 



u 0.48 



0.46 



0.44 



_i 1 1 1 i_ 



n, =1.9-10 




best fit: 
e 4 =-0.532(2)-KT* 



3 4 



5 6 7 8 



0.48 



0.44- 



_i I i_ 



best fit: 

e =-0.38(1)-10-' 




(a) 

Fig. 2. (a) Same as Fig. 1 but for the 4 th cumulant. (b) Data for the derivative of the 
5 th cumulant. Open triangles denote results (not used in the fit) which have a \ 2 > 1-5 
indicating that successive Monte-Carlo iterations are not consistent with each other. 



The best fit values for a and e$ displayed in Figs. 2 a, b are still preliminary 
as a more detailed error analysis has to be made. Also for the n = 5 case 
the Monte-Carlo statistics should be improved. Note that each high-statistic 
point in Fig. 2 b took about 30 days runtime on a Xeon 3.0 GHz machine. 
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4. Summary and Outlook 

• Two additional pcrturbative coefficients e^, e$ for the polaron g.s. 
energy have been determined by the method of Bogolyubov and 
Plechkov (rediscovered independently). This amounts to perform- 
ing a 4-loop and 5-loop calculation in Quantum Field Theory. 

• The method is based on a combination of Monte-Carlo integration 
techniques and deterministic quadrature rules for finite (5 (tem- 
perature) and on a judicious extrapolation to (3 — > oo (zero tem- 
perature). As a check the value of calculated analytically by 
Smondyrev has been reproduced with high accuracy. 

• The cancellation in n th order is not among many individual dia- 
grams but among the much fewer terms in the integrand of the 
(3n — 3)-dimensional integral (see Eq. (9)). 

• The method can be simply extended to the calculation of higher- 
order terms in the small-coupling expansion of the effective mass 
m*(a) for a moving polaron. 

• Generalizing this approach to relativistic QFT in the worldline 

1 9 

representation 1 ^ and calculation of higher-order terms for the 
anomalous magnetic moment of the electron is under investiga- 
tion. New challenges arise from the divergences which now occur 
and the need for renormalization. 
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